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' Abstract. In this second paper in a series dedicated to developing efficient numerical techniques for the deblurring 

, Cosmic Microwave Background (CMB) maps, we consider the case of asymmetric point spread functions (PSF). 

' Although conceptually this problem is not different from the symmetric case, there are important differences from 

' the computational point of view because it is no longer possible to use some of the efficient numerical techniques 

that work with symmetric PSFs. We present procedures that permit the use of efficient techniques even when this 



Q ' condition is not met. In particular, two methods are considered: a procedure based on a Kronecker approximation 

technique that can be implemented with the numerical methods used with symmetric PSFs but that has the 
C/3 ■ limitation of requiring only mildly asymmetric PSFs. The second is a variant of the classic Tikhonov technique 



i-i 



X 



that works even with very asymmetric PSFs but that requires discarding the edges of the maps. We provide 
details for efficient implementations of the algorithms. Their performance is tested on simulated CMB maps. 

Key words. Methods: data analysis - Methods: statistical - Cosmology: cosmic microwave background 

1. Introduction sources will vary over the sky. As a result, even a suc- 

I 1 1 1 cessful separation of the components contributing to the 

In the first paper in the series | Vio et al. | mm (VNT) j^icrowave radiation will provide results of inhomogeneos 

have stressed the advantage of deblurrmg small patches of Therefore, even if some characteristics of CMB are 

the sky m cosmic microwave background (CMB) studies: estimated on full sky maps, it will be important to check 

first, It helps to recover high frequencies smoothed out by ^hese results on smaller sky patches where CMB largely 

the instrument's PSF. Second, a better understanding of dominates the other components (i.e., no component sep- 

sky emissions, from foregrounds m particular, is achieved ^^^^ion is necessary as, e.g., at high Galactic latitude and 

if multifrequency sky maps are compared on a common high observing frequency) and data are free from in- 

resolution. Third, some map-based component separa- gtrumental and/or observational problems. We stress that 

tion algorithms, such as independent component analysis ^ ^^-^ resolution is always possible also when maps of 

jBaccigalupi et al. || 2000t | Maino et al. || 2002i) , require in- ^^g^,^^^ resolutions are composed together, for example, 

put maps with similar level of degradation. Furthermore, ^^im^g^ the method described in iTegmarki il999lb 
although the aim of satellite missions such as Planck and 

MAP is to obtain full sky maps of the CMB, the strength VNT suggested a deblurring approach based on 

of the CMB over other backgrounds or contaminating Tikhonov regularization whose computational cost is com- 

parable to that of classic frequency domain methods but 

Send offprint requests to: R. Vio that leads to more reliable and stable deblurring estimates. 
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These new efficient implementations make the Tikhonov 
technique a promising tool for deblurring CMB maps. The 
main limitation of the method is the requirement of sym- 
metric PSFs, which prevents its application in the general 
experimental context with no restriction on the form of 
the PSF. This point is especially important in CMB ex- 
periments where an asymmetric PSF can alter the results 
of the analyses, for example, by distorting the estimated 
angular power spectrum or by altering the measure of the 
degree of nonCaussianity of the maps. Since some exper- 
imental situations are not well controlled, one should de- 
velop algorithms that can cope with the worst possible 
scenario; it is thus important to develop deblurring meth- 
ods for asymmetric PSFs. 

In Sect. 121 we formalize the problem and propose two 
solutions in Sects. 01 and ^ the first one is a very effi- 
cient method but requires only mildly asymmetric PSFs, 
whereas the second, although less efficient, is not limited 
by the particular form of the observing beam. A modifi- 
cation of the second method is considered in Sect.[5l The 
results of numerical simulations to test the performance of 
the different methods are presented in Sect. El In Sect.H 
we close with final comments and conclusions. 



2. Formalization of the problem 

We make use of the same formalism adopted in VNT. 
When a two-dimensional object /(f , 77) is observed 
through an optical invariant (linear) system, it is seen as 
an image g{x,y), 



hix~^,y-i])f{^,r]) dr], (1) 



where the space-invariant point-spread function (PSF) 
h{x — (,,y — Tj) represents the blurring action of the op- 
tical instrument. This model is only theoretical, in practi- 
cal applications we only have discrete noisy observations 
of the image g{x, y), which we model as a discrete linear 
system 

g = Hf + z, (2) 

where: g = vec(G) and / = vec(i^) are one-dimensional 
column arrays containing, respectively, the observed image 
G and the true images F in stacked order, z is an array 
containing the noise contribution (assumed to be addi- 
tive), and H is a matrix that represents the discretized 
blurring operator. 

There are two problems in obtaining an estimate of / 
from g: the size of the matrix H, which is large even for 
moderate size images, and the very ill-posed nature of the 
problem. VNT proposed a deblurring method for CMB 
applications that is efficient and numerically stable; it is a 
Tikhonov regularization approach where the estimate 
of / is defined as 



with A a scalar regularization parameter. L is often the 
identity matrix or a discrete derivative operator of some 
order. 

An additional problem is the selection of boundary 
conditions (BC) to account for data outside the field of 
view. VNT found that better deblurring estimates of CMB 
maps were obtained with reflexive BC. This choice leads 
to reliable and stable regularization parameter estimates, 
and helps suppress spurious features such as Gibbs oscil- 
lations. Periodic and zero BC impose edge discontinuities 
which bias the image Fourier coefficients and affect the 
reliability of the regularization parameter estimates ob- 
tained through generalized cross-validation (GCV). 

On the other hand, efficient implementations of 
Tikhonov deblurring with reflexive BC require symmetric 
PSF (although not necessarily separable), i.e., h{x,y) = 
h{—x, y) — h{x, —y) ~ h{~x, —y). Since the PSF may be 
asymmetric in some practical applications (there are indi- 
cations that this may be the case for PLANCICs optics), 
we consider efficient implementations of Tikhonov deblur- 
ring that can be used in this case. In particular, in Sect.O 
we consider a very efficient method based on reflexive BC 
that, however, works only when the PSF is slightly asym- 
metric. In Sect. 01a less efficient method, based on periodic 
BC, is presented; its performance is not sensitive to the 
specific form of the PSF. For this last case we have to 
modify the traditional GCV to provide regularization pa- 
rameter estimates in the presence of edge discontinuities. 

3. Kronecker approximation 

One possible alternative when the PSF is not symmetric 
is to determine if the PSF is at least separable; that is. 



h{x,y) = K(x)hy(y) . 



(4) 



If this is the case, then the corresponding x matrix 
H can be written as 



H ^ A®B. 



(5) 



where A and B are nxn matrices, and (g) is used to denote 
a Kronecker product 



A®B = 



( aiiB ai2B 
a2iB 022-6 

\ OniB a„2-B 



ainB ^ 
a2nB 



.B j 



argmin(||i?/-g||2 + A2||i/ 



(3) 



For such a structured matrix, algorithms can be imple- 
mented efficiently. The cost is O(n^), which is slightly 
more than the O(n^logn) required of transform-based 
methods, but is still very reasonable for large images; see 
VNT for further details. Therefore, we should exploit this 
structure if we can recognize that the blur is separable. 

The key point that permits the development of an ef- 
ficient algorithm is that if P is an n x n array of pixels 
representing the observed PSF and the blur is separable. 
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then P is a rank-one matrix; that is, P can be written as 
an outer product 



P = ab' 



(6) 



where a and b are n x 1 vectors. For reflexive boundary 
conditions the matrices A and B (where P is such that 
Pij is the center of the PSF) are Toephtz-plus-Hankel of 
the form 



/a, ■■ ■ ai 



A = 



ai 



V a-a ■■■ ai ) \ ax ■■ ■ a^^i / 

(7) 



and 



B = 



V 



bn 



hi 



•j+i ■ ■ ■ 



bn \ 



bi 



V 



bi 



(8) 

If the PSF if effectively separable, then a and b can be 
obtained, respectively, by discretizing the functions hx{x) 
and hy{y) in Eq. (gj). 

Of course, these arguments do not hold when the PSF 
is not separable since it cannot be expressed exactly as an 
outer product of two vectors. However, if the PSF is only 
approximately non-separable, then it is possible to work 
with a separable approximation of P. The most natural 
approach is to make use of the singular value decomposi- 
tion (SVD) of P: 



(9) 



where the singular values Ui satisfy tii > • • • > (t^ > 0, r is 
the rank of P, and the singular vectors Ui and Vi satisfy 



T 



T 



1 if i = j 
if*^j 



(10) 



The best rank-one approx imation of P is given by (e.g. 
IColub Va,n T.oari iFmfl^) 



(11) 



Thus, we may take a 
implies 

P « ab^ and 



y^Ui and b = ^/oiVi, which 
HKA(g>B. (12) 



Although this often works well in practice, it is generally 
preferable, and mathematically more satisfying, to have 
an optimal approximation. Specifically, we would like to 
find a Kronecker product A® B that minimizes 



minlliJ- A(g)B| 



(13) 



where 1 1 • 1 1 is a chosen norm, and the minimization is 
done over all Kronecker products A^ B. Approximations 
in the case of the Frobenius no rm ( || • || f) have re- 
cently received a lot of attention. IVan Loan fc Pitsianis I 
(199i|) developed the idea for general matrix approxima- 
tions, which was made computationally efficient for im- 
ag e restoration problems with zero boundary conditions 
bv lKamm fc Naevl l)200(]|) . For refiexive boundary condi- 
tions, the optimal approximation can be computed using 
the follow ing theorem, which was recently established by 
iNaev et a l. ( 2002). 

Theorem 1. Let P be an nxn PSF. For reflexive bound- 
ary conditions: 



\H-A®B\ 



ab \\f 



where P = RPR , d = Ra, b = Rb and R is the 

Cholesky factor of the nxn symmetric Toeplitz matrix 
with first row [nlOlOlOl---]. 

The proof of this theorem, which requires many tedious 
details, is given in Nag y _et a l. (2002). Based on this theo- 
rem, an algorithm for constructing the optimal Kronecker 
product approximation is as follows: 
Algorithm: To construct the approximation H « A^B: 



Compute R 
Construct Pr = RPR^ 
Compute the SVD: Pr 
Construct the vectors: 



J2 cfkUkvl 



a — R ^1*1 and b = \fa\ R 

— Construct the matrices A and B from a and b (as 
described above). 

If the PSF image array is of size m x m, then the cost 
of constructing this optimal Kronecker product approx- 
imation is only O(to^), which is relatively cheap if the 
width of the PSF is small compared to the dimension of 
the blurred image (i.e., m <^ n). We use this scheme in 
our computations because it produces a provably optimal 
approximation, and because it is essentially equivalent in 
cost to the straight forward approach given in Eq. (|12|l . 

4. Image windowing 

The Kronecker approximation method performs poorly 
when the PSF is very asymmetric. We now consider a 
deblurring method based on periodic BC which can be 
implemented using the fast Fourier transform. 

A problem with periodic BC is the effect of edge dis- 
continuities in regularization parameter estimates. This 
effect can be reduced by considering an average GCV de- 
fined by an estimated spectrum. More precisely, the GCV 
in the spectral domain is (see VNT) 



GCV(A) 



^2 + 



(14) 
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where: {ai} and {Si} are the eigenvalues of H and L, re- 
spectively, g — ^g, with ^ the bidimensional Fourier ma- 
trix, and n is the number of pixels in the image. To reduce 
edge effects on g we use an estimate of || g |p (the spec- 
trum) based on bidimensional windowing of the Fourier 
transform: g is estimated via gi = ^{wi Q g), where " 
© " denotes element- wise matrix multiplication and Wi 
is a bidimensional function (window) tapered smoothly to 
zero at the image edges. A typical choice for wi is the 
Hanning window (see Fig.^l, but many other alternatives 
are available. The need for windowing is evident in Fig. 12 
which shows a simple example based on the realization 
of a pure sinusoidal process. It is clear that without win- 
dowing the (classic) estimated power-spectrum of the sig- 
nal strongly depends on the the specific sampling pattern. 
Note the spread of power at all the frequencies, visible in 
panels (a) and (d), which has deleterious effects in regu- 
larization parameter estimates. Windowing stabilizes the 
power-spectrum estimates. 

Once the parameter A is estimated, a second window- 
ing is necessary to also reduce edge effects in the deblur- 
ring stage, i.e., the deblurring operation has to be carried 
out on a windowed map ^2 — ^{^2 Q g)- In this step, 
however, W2 should distort the map as little as possible. 
A possible solution is a window that does not alter the 
image within a central subimage. We use the following 
modification of the classic Banning window 



0.25 X ax (3 



1 < hj < N^; 

N-N^<i,j<N; (15) 
otherwise; 



where a = [1 — cos(7r {i — l)/iV^)] , and /3 ~ [1 — cos(7r {j — 
1)/A^u))] so that the pixels in the central subimage are not 
modified and the image has continuous first derivatives at 
the edges (see Fig. QJ. This window approaches the clas- 
sic two-dimensional Hanning window as N/2 and 
tends to the rectangular windows as A^^i 0. The pa- 
rameter iV^, thus determines the filtering characteristics, 
in particular the frequency pass-band, of the window. Its 
"optimal" value depends on many factors such as the noise 
level, the form of the PSF and the specific realization of 
the process. 

5. Reflexive image extension 

Another method that allows, at least in principle, the use 
of reflexive BC with asymmetric PSFs and which does 
not require discarding any data, consists of extending the 
original image X according to the scheme 



X X, 

Xr Xrr 



where X^ is obtained by "flipping" the columns of X, Xr is 
obtained by "flipping" the rows of X, and Xrc is obtained 
by "flipping" the rows and columns of X. Since the result- 
ing image is periodic, periodic BC can be used without in- 
troducing discontinuities. However the computational cost 




10 20 30 40 50 60 70 80 90 100 

Fig. 1. Slice comparison of the two-dimensional classi- 
cal Hanning and rectangular windows with the modified 
Hanning window {N^ — 20). 



a) N=191, f =0.15 Hz 
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Fig. 2. Power-Spectra of four different realizations of the 
process x[t] = sin(27r/ot), t = 0, 1, . . . , TV for /o = 0.150 
and 0.157 Hz, N = 191 and 200. The true power-spectrum 
(not shown) is a 5- function centered at /q. 



is higher as the new image is nine times larger. A more 
feasible approach is to make Xrc, Xc, and Xr only part of 
the image, say bands of thickness A^cxt of the same order as 
the width of the PSF. Since the discontinuities introduced 
by the BC are expected mainly at the edges, this image 
extension keeps edge effects on a part of the image that is 
later discarded. This method, however, has the disadvan- 
tage of not providing full reflexive BC and thus it has to 
be implemented using a windowing operation similar to 
that presented in Sect ^ with iV^, = iVcxt- 
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S/N 


Wiener 


Ref. Ext. 


Tikhonov 


Wind. Ref. Ext. 


Wind. Tikh. 


Kron. Approx. 


rrms (%) 


rrms {%)^ 


rrms (%) 


rrms (Vo)^ 


rrms (%) 


rrms (%)^ 


2 


50.71 ±0.12 


51.09 ±0.20 


50.18 ±0.21 


50.99 ±0.18 


49.99 ± 0.20 


52.37 ±0.37 


10 


45.64 ± 0.07 


52.01 ±0.39 


47.08 ± 0.28 


50.90 ± 0.37 


45.50 ±0.13 


* 


100 


41.29 ±0.04 


* 


* 


* 


42.83 ±0.19 


* 



Table 1. Summary of the results obtained for a simulated sky map contaminated with 100 different realization of a 
white noise process at the PLANK-LFl frequencies and axial ratio of the elliptical PSF equal to 1 : 1.3 (see text). The 
central 300 x 300 pixels of the images are considered. The original maps consisted of 364 x 364 pixels, i.e., a border 
of 32 pixels has been removed from each side of the images, corresponding to about four times the dispersion of the 
PSF along the major axis. — 32 for the modified Hanning window used for the windowed reflexive expansion and 
the windowed Tikhonov methods and for the calculation of the regularization parameter in the reflexive expansion 
and Tikhonov methods (see text). All of the methods have adopted a discrete Laplacian for L and ( except for the 
Kronecker approximation) periodic BCs. The Kronecker approximation uses reflexive BCs. The relative root mean 
square (rrms) is defined as the ratio of the residual root mean square (rms) to the rms of the true signal. The asterisk 
means unstable results. 

Calculated on full 364 x 364 pixel images. 



S/N 


Wiener 


Ref. Ext. 


Tikhonov 


Wind. Ref. Ext. 


Wind. Tikh. 


Kron. Approx. 


rrms (%) 


rrms (%)^ 


rrms (%) 


rrms {%)^ 


rrms (%) 


rrms (%)^ 


2 


48.42 ± 0.12 


49.21 ±0.18 


47.98 ±0.17 


48.97 ±0.16 


47.91 ±0.17 


56.82 ±0.72 


10 


43.35 ± 0.07 


52.50 ±0.50 


44.56 ±0.22 


50.86 ± 0.45 


43.19 ±0.13 


* 


100 


38.90 ± 0.04 


* 


* 


* 


39.67 ±0.14 





Table 2. As in Table with the axial ratio of the elliptical PSF equal to to 1 : 2 



S/N 


Wiener 


Ref. Ext. 


Tikhonov 


Wind. Ref. Ext. 


Wind. Tikh. 


Kron. Approx. 


rrms (%) 


rrms {%)^ 


rrms (%) 


rrms (Vo)^ 


rrms (%) 


rrms (%)^ 


2 


51.04 ±0.55 


51.91 ±0.56 


50.64 ±0.62 


51.85 ±0.55 


50.60 ± 0.63 


51.17 ±0.65 


10 


45.73 ±0.38 


52.81 ± 1.15 


47.11 ±0.74 


51.91 ± 1.07 


45.56 ±0.41 


* 


100 


41.12 ±0.26 


* 


* 


* 


42.65 ± 0.38 


* 



Table 3. As in Table with the only difference that a new Gaussian random field is generated for each simulation. 



See Attached Fig. 3 

Fig. 3. Grayscale image of the central 290 x 290 pixels 
of a simulated sky map at the PLANCK-LFI frequencies 
with S/N = 2 and axial ratio of the elliptical PSF equal to 
1 : 1.3 (see text). The original map was of 354 x 354 pixels, 
i.e., a border of 32 pixels has been removed from each 
side of the image, corresponding to about four times the 
dispersion of the PSF along the major axis. iV^, = 32 has 
been used for the dewindowed modified Hanning methods 
as well as Periodic BC and discrete Laplacian for L. For 
the approximated Kronecker method a reflexive BC has 
been adopted. 



See Attached Fig. 4 

Fig. 4. As in Fig. 01 but with S/N = 10 and axial ratio 
= 1 : 1.3. 

See Attached Fig. 5 

Fig. 5. As in Fig. |21 but with S/N = 2 and axial ratio 
= 1:2. 

6. Numerical experiments 

6.1. A preliminary experiment 

To support our arguments and study the performance of 
the approximations and methods described in the previ- 
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S/N 


Wiener 


Ref. Ext. 


Tikhonov 


Wind. Ref. Ext. 


Wind. Tikh. 


Kron. Approx. 


rrms (%) 


rrms {%)^ 


rrms (%) 


rrms (Vo)^ 


rrms (%) 


rrms (%)^ 


2 


48.63 ± 0.50 


49.99 ±0.51 


48.32 ±0.57 


49.75 ± 0.50 


48.28 ± 0.57 


57.20 ± 1.11 


10 


43.23 ± 0.33 


52.61 ± 1.32 


44.12 ±0.54 


51.52 ± 1.13 


43.10 ±0.36 


* 


100 


38.62 ±0.22 


* 


* 


* 


39.33 ±0.27 


* 



Table 4. As in Table [21 with the only difference that a new Gaussian random field is generated for each simulation. 



Windowed Tikhonov 



See Attached Fig. 6 

Fig. 6. As in Fig. |3|but with S/N = 10 and axial ratio 
= 1:2. 



ous sections, we present the result of numerical simulations 
based on a Gaussian random process with statistical char- 
acteristics similar to those expected for the CMB emission 
at the frequencies typical of PLANCK-LFl (see VNT). 
The PSF is Gaussian with elliptical symmetry, the FWHM 
along the major axis is « 8 pixels (about two times the 
worst spatial resolution expected for PLANCK), and with 
axes forming an angle of 45° with the edges of the map. 
We consider the S/N ratios 2, 10, 100; and two values of 
the axial ratio: 1 : 1.3 and 1 : 2. The axial ratios and the 
width of the PSF used in the experiment are by far less 
favorable than those expected for PLANCK and repre- 
sent a sort of worst possible scenario. The reason for using 
S/N ratios much more larger than the value expected for 
PLANCK (« 2) is that some characteristics of the method 
are hidden by high noise contamination. Six deblurring 
methods are tested: classic Wiener, Tikhonov (with peri- 
odic BC and discrete Laplacian for L) applied both to g 
and g (windowed Tikhonov), reflexive extension method 
(with periodic BC and discrete Laplacian for L) applied 
both to a windowed and to an unwindowed extended im- 
age, and the Kronecker approximation with reflexive BC. 
Since the random field is Gaussian and stationary and the 
noise is assumed white, classical Wiener filtering is ex- 
pected to provide the smallest mean square error among 
linear filters. However, this filter requires knowledge of 
the spectrum of the unknown signal which is not available 
in practice. We use Wiener deblurring, based on the real 
spectrum of the signal, as a sort of benchmark to assess 
the performance of Tikhonov methods. For this reason, 
it has also been implemented in a way that avoids edge 
effects. 

The simulations have been conducted under two dif- 
ferent scenarios. We first fix the sky and generate different 
realizations of the noise process. Then, to account for the 
variability of the random field, we simulate different real- 
izations of the random field and the noise process. 

Tabs 1 114] and Figs l3l6l confirm that the Kronecker ap- 
proximation method docs not perform well with a very 
asymmetric PSF. Furthermore, it is evident that the win- 



10 20 30 40 50 



10 20 30 40 50 



Fig. 7. Standard deviation (7^ of the residual between the 
deblurred and the true maps vs. the width Nb of the re- 
moved border for the map in Fig. El 



dowed Tikhonov method performs the best, close to the 
Wiener filter, even in the case of very high S/N that is 
troublesome for the other methods. Figs. I7I9I show the 
typical behavior of the standard deviation of the resid- 
uals of the deblurred and the true maps as larger bor- 
ders are removed from the frames. The effect of the Gibbs 
phenomenon is evident in the figures, especially those ob- 
tained by directly deblurring the unwindowed g and for 
high values of S/N. These figures indicate that for mod- 
erate S/N ratios a border of thickness 3-4 times the dis- 
persion of the PSF has to be removed after the deblurring 
to reduce edge effects. By increasing N^^, this method still 
performs better even for very high S/N ratios. In typical 
CMB applications the S /N is low, thus here we do not con- 
sider the question of finding an "optimal" value of A^u, . In 
fact, our simulations indicate that a value of equal to 
3-4 times the dispersion of the PSF along the major axis is 
a reasonable choice. This value corresponds approximately 
to the thickness of the border in which the blurred image 
is influenced by data outside of the field of view. 

Finally, as shown in Tables 1 1141 the methods provide 
similar results for low S/N ratios, especially when a suf- 
ficiently large number of edge pixels is removed from the 
image. This is not surprising as edge effects remain close 
to the edges and high noise levels hide the effects of the 
smallest eigenvalues of H that cause the ill-posedness of 
the deblurring operation. 

Tables|Sland|H|show the results of 100 simulations simi- 
lar to those presented previously but under the conditions 
expected for PLANCK-LFl. As expected, given the low 
S/N ratio, the performance of the various methods is sim- 
ilar. 
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FWHM (arc.min.) 


Wiener 


Windowed Tikhonov 


Kronecker Approx. 


rrms (%) 


rrms (%) 


A 


rrms (Vo)^ 


A 


10 


30.36 ± 0.06 


30.35 ±0.06 


0.69 ±0.02 


30.45 ± 0.05 


0.69 ± 0.01 


14 


32.64 ± 0.07 


36.68 ±0.09 


0.79 ± 0.04 


37.96 ± 0.08 


0.79 ± 0.04 


23 


36.74 ± 0.08 


36.68 ± 0.09 


0.79 ± 0.04 


36.96 ± 0.08 


0.79 ± 0.02 


33 


40.54 ±0.10 


40.47 ±0.11 


0.77 ±0.04 


40.70 ± 0.11 


0.79 ± 0.02 



Table 5. Summary of the results concerning the deblurring of a sky map contaminated with 100 realizations of a 
Gaussian random process whose statistical properties are similar to those expected of the CMB sky observed with 
four channels of PLANCK-LFl for beams with elliptical symmetry. The central 320 x 320 pixels of the images are 
considered. The original maps consisted of 350 x 350 pixels (corresponding to a sky area of about 20° x 20°), i.e., a 
border of 15 pixels has been removed from each side of the images. The techniques used are Wiener filtering, windowed 
Tikhonov (with periodic boundary conditions and discrete Laplacian for L) and the Kronecker approximation (with 
reflexive boundary conditions and discrete Laplacian for L). Here S/N — 2, FWHM is the full width at half maximum 
along the major axis of the PSF, the axial ratio for the PSF is 1 : 1.3 and iV^, = 15. The relative root mean square 
(rrms) is defined as the ratio of the residual root mean square (rms) to the rms of the true signal. When applicable, 
the mean values and dispersions of the GCV estimates of A are also shown. 

^ Calculated on the entire images of 350 x 350 pixels. 



FWHM (arc.min.) 


Wiener 


Windowed Tikhonov 


Kronecker Approx. 


rrms (%) 


rrms (%) 


A 


rrms (%)'' 


A 


10 


30.68 ±0.26 


30.71 ±0.29 


0.73 ± 0.03 


30.73 ± 0.28 


0.73 ± 0.03 


14 


32.68 ±0.24 


32.68 ±0.28 


0.75 ± 0.04 


32.73 ±0.27 


0.76 ± 0.04 


23 


36.74 ± 0.24 


36.67 ±0.28 


0.78 ± 0.04 


36.78 ± 0.26 


0.79 ± 0.04 


33 


40.65 ±0.27 


40.52 ±0.30 


0.77 ±0.05 


40.73 ± 0.29 


0.81 ±0.05 



Table 6. As in Table El with the only difference that a new Gaussian random field is generated for each simulation. 



Windowed Tikhonov 



Windowed Til^tionov 




10 20 30 40 50 



Fig. 8. Standard deviation ar of the residual between the 
deblurred and the true maps vs. the width Nb of the re- 
moved border for the map in Fig. 



6.2. Numerical simulations with realistic CMB maps 

The next step is to test the performance of the deblurring 
procedure with realistic CMB maps from the point of view 
of the angular power spectrum. As usual, it is defined as 
the set of coefficients Ce of the two-point correlation func- 
tion expanded in Legendre polynomials; £ marks the power 
at the angular scale given approximately by ~ 180/^. We 
compare the quality of the reconstruction with the corre- 
sponding one in VNT, which assumed a circular Gaussian 
beam. 



Fig. 9. As in Figs. and |H1 but with S/N = 100 



We consider the same templates as in VNT: the re- 
gion is a squared patch (350 x 350 pixels) with side of 
about 20°, centered at I = 90° , b = 45° (Galactic co- 
ordinates). The latitude is high enough that CMB emis- 
sion dominates over for egrounds, assumed to be repre- 
s ented by synchrotr on l)Haslam et al. 1 1 198^ and dust 
fSch leeel et al. IITqQSI) emission. We neglect contributions 
of point sources. The CMB model, in agreement with 
current experimental results Jd e Berna rdis et al. I l2002t 
iHalverson et al. Il2002t iLee et al. ■ ■2001). corresponds to 
a flat Friedmann-Robertson- Walker (FRW) metric with 
a cosmological constant (70% of the critical density), 
Hubble parameter today Ho = lOOh km/sec/Mpc with 
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h = 0.7 baryons at 5% and Cold Dark Matter (25% CDM), 
with a scale-invariant Gaussian initial spectrum of adia- 
batic density perturbations. 

The PLANCK-LFl instrument works at frequencies 30, 
44, 70, and 100 GHz. We assume nominal noise and angu- 
lar resolution corresponding to the four frequencies 30, 44, 
70, and 100 GHz at which the PLANCK-LFl instrument 
works. The simulated maps are blurred through Gaussian 
PSF's with elliptical symmetry. In particular, the follow- 
ing FWHM's along the major axis have been used: ~ 33' 
at 30 GHz, « 23' at 44 GHz, « 14' at 70 GHz, « 10' at 
100 GHz, with axial ratio set to 1 : 1.3, and axes forming 
an angle of 45° with the edges of the map. Simulated white 
noise, with rms level as expected for the considered chan- 
nels, has been added to the maps. Since we choose to work 
with a pixel size of about 3.5 arcminutes, the noise rms 
are .042, .049, .042 and .043 mK in antenna temperature 
at 30, 44, 70, 100 GHz, respectively. 

Again, as in VNT, we can see two important char- 
acteristics of the deblurring process: first it reconstructs 
the correct shape and amplitude of the part of the spec- 
trum which is mildly affected by the PSF; second, it re- 
constructs part of the power where the signal is degraded 
substantially by the PSF and noise. In fact, we see that 
the performance of the deblurring method is almost the 
same as that reported in VNT for symmetric PSF. Indeed, 
in the 30, 44, 70 and 100 GHz cases the spectrum is re- 
constructed up to £ ~ 400, 500, 700 and 800 respectively. 

We conclude that, for what concerns the power spec- 
trum analysis, the proposed technique seems to work well 
also with asymmetric shapes of the instrumental beam. 

Note also that for a fixed A the Tikhonov estimate Q 
is a linear function of the data and therefore the covariance 
matrix of the deblurrcd field can be obtained by propagat- 
ing the covariance matrix of the original field through the 
linear operators. Even if the GCV estimate of A makes the 
estimate nonlinear, for large samples the linear approxi- 
mation that assumes A is fixed is reasonable. In any case, it 
is necessary to stress that the most feasible way to account 
for the changes in the statistics of the CMB through any 
estimation process (to include the effects of, for example, 
edge effects, possible unremoved point sources, asymme- 
try and/or non-stationary PSF etc.), is probably through 
Monte Carlo simulations where one applies the same al- 
gorithm to the data and to simulated maps to compare 
and fit the best cosmological model. The low computa- 
tional cost of our algorithm is certainly interesting with 
this respect. 

7. Discussion and conclusions 

We have considered Tikhonov regularization for deblur- 
ring CMB maps in real space. As shown in VNT, this 
approach permits the development of algorithms that are 
more flexible and robust than those based on frequency- 
space methods. The methods developed in VNT, however, 
apply to the case of symmetric PSFs for which efficient 
methods can be implemented with reflexive BC. In the 



original blurred (eiiiptical baam) 
i 0.007 t I 




Fig. 10. Angular power spectrum at 30 GHz in different 
steps of the analysis. 

present paper we have considered the more general case 
of asymmetric PSFs. 

We have presented a method based on a Kronecker sep- 
arable approximation of the PSF that can be used with 
mildly asymmetric PSFs. For more asymmetric cases we 
presented a periodic BC approach that can be efficiently 
implemented with image windows and fast Fourier trans- 
forms. Windowing is necessary to reduce edge effect in the 
selection of the regularization parameter. Of course, one 
can easily derive a GCV function that only takes into ac- 
count pixels away from the edges but the computational 
cost is higher. 

We have applied our methodology to simulated skies 
at typical CMB frequencies. We considered test signals 
with known statistics, as well as realistic simulations of 
the CMB sky contaminated by noise whose rms is that 
expected for the low frequency instrument aboard the 
PL^A^Cif satellite. This case is particularly interesting for 
application of a deblurring procedure, as the instrument 
observes the sky at 30, 44, 70 and 100 GHz with very dif- 
ferent PSFs of resolution 33, 22, 14, 10 arcminutes. We 
found that the proposed methodology performs as well or 
better than the Wiener benchmark that relies on the true 
spectrum and that avoids edge effects. 
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Fig. 11. Angular power spectrum at 44 GHz in different Fig. 13. Angular power spectrum at 100 GHz in different 
steps of the analysis. steps of the analysis. 
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Fig. 12. Angular power spectrum at 70 GHz in different 
steps of the analysis. 
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